/**
 * @file EdgeDetection.cpp
 * @author enemy1205 (enemy1205@qq.com)
 * @brief Canny边缘检测算法实现
 * @date 2021-09-27
 */
#include"EdgeDetection.h"
#include <iostream>

//阶乘
//int factorial(int n){
//    int fac = 1;
//    //0的阶乘
//    if (n == 0)
//        return fac;
//    for (int i = 1; i <= n; ++i){
//        fac *= i;
//    }
//    return fac;
//}
//阶乘递归写法
float factorial(int n) {
    if (n == 0)
        return 1;
    else return static_cast<float>(n) * factorial(n - 1);
}

/**
 * @brief 计算平滑算子
 * @note N!/(K!*(N-K)!)
 * @example 5阶：[1, 4, 6, 4, 1]
 */
Mat _Sobel::calSmoothOperator(int k_size) {
    Mat smooth_operator = Mat::zeros(Size(k_size, 1), CV_32FC1);
    for (int i = 0; i < k_size; ++i) {
        smooth_operator.ptr<float>(0)[i] = factorial(k_size - 1) / (factorial(i) * factorial(k_size - i - 1));
    }
    return smooth_operator;
}

/**
 * @brief 计算差分算子
 * @example [1, 2, 0, -2, -1]
 */
Mat _Sobel::calDifferOperator(int k_size) {
    Mat differ_operator = Mat::zeros(Size(k_size, 1), CV_32FC1);
    //差分算子的前身是n-2阶的二项式展开式，与平滑算子相差1
    Mat pascalSmooth_previous = calSmoothOperator(k_size - 1);
    for (int i = 0; i < k_size; ++i) {
        if (i == 0) differ_operator.ptr<float>(0)[i] = 1;//最左列置1
        if (i == k_size - 1) differ_operator.ptr<float>(0)[i] = -1;//最右列置-1
        else
            differ_operator.ptr<float>(0)[i] =
                    pascalSmooth_previous.ptr<float>(0)[i] -
                    pascalSmooth_previous.ptr<float>(0)[i - 1];//其余列=当前列平滑算子-前一列
    }
    return differ_operator;
}

/**
 * @brief 二维图像卷积
 * @param src 源图
 * @param dst 结果图
 * @param kernel 卷积核
 * @param ddepth 深度
 * @note ddepth 似乎仅CV_32FC1可行
 */
void _Sobel::conv2D(Mat &src, Mat &dst, Mat &kernel, int ddepth) {
    Mat kernelFlip;//水平垂直翻转
    flip(kernel, kernelFlip, -1);//水平垂直翻转
    filter2D(src, dst, ddepth, kernelFlip, Point(-1, -1), 0, BORDER_DEFAULT);
}

/**
 * @brief 可分离卷积———先垂直方向卷积，后水平方向卷积
 * @param src 源图
 * @param dst 结果图
 * @param kernel_Y 垂直方向卷积核
 * @param kernel_X 水平方向卷积核
 * @param ddepth 深度
 */
void _Sobel::sepConv2D_Y_X(Mat &src, Mat &dst, Mat &kernel_Y, Mat &kernel_X, int ddepth) {
    Mat dst_kernel_Y;
    conv2D(src, dst_kernel_Y, kernel_Y, ddepth); //垂直方向卷积
    conv2D(dst_kernel_Y, dst, kernel_X, ddepth); //水平方向卷积
}

/**
 * @brief 可分离卷积———先水平方向卷积，后垂直方向卷积
 * @param src 源图
 * @param dst 结果图
 * @param kernel_X 水平方向卷积核
 * @param kernel_Y 垂直方向卷积核
 * @param ddepth 深度
 */
void _Sobel::sepConv2D_X_Y(Mat &src, Mat &dst, Mat &kernel_X, Mat &kernel_Y, int ddepth) {
    Mat dst_kernel_X;
    conv2D(src, dst_kernel_X, kernel_X, ddepth); //水平方向卷积
    conv2D(dst_kernel_X, dst, kernel_Y, ddepth); //垂直方向卷积
}

/**
 * @brief 实现opencv中sobel效果，获得梯度图像
 * @param src 待处理灰度图
 * @param k 算子尺寸
 */
_Sobel::_Sobel(Mat &src, int k) {
    L2graydient = false;
    Mat smooth_operator = calSmoothOperator(k);//平滑算子
    Mat differ_operator = calDifferOperator(k);//差分算子
    Mat smooth_operator_t = smooth_operator.t();//转置处理
    Mat differ_operator_t = differ_operator.t();
    //可分离卷积———先垂直方向平滑，后水平方向差分——得到垂直边缘
    sepConv2D_Y_X(src, dst_X, smooth_operator_t, differ_operator, CV_32FC1);
    //可分离卷积———先水平方向平滑，后垂直方向差分——得到水平边缘
    sepConv2D_X_Y(src, dst_Y, smooth_operator, differ_operator_t, CV_32FC1);
    //计算梯度幅值:精确计算or近似计算
    if (L2graydient)
        magnitude(dst_X, dst_Y, edgeMag); //开平方
    else
        edgeMag = abs(dst_X) + abs(dst_Y); //绝对值之和近似
    //边缘强度（近似）
    edgeMag = abs(dst_X) + abs(dst_Y);
    convertScaleAbs(edgeMag, edgeMag); //求绝对值并转为无符号8位图
}


/**
 * @brief 构造函数
 * @param path 源图路径
 * @param kernel_size 算子尺寸
 */
EdgeDetection::EdgeDetection(const std::string &path, int kernel_size) {
    this->src = imread(path);
    this->k_size = kernel_size;
//    padded_img = Mat::zeros(Size(src.cols + 2 * (k_size / 2), src.rows + 2 * (k_size / 2)), CV_8UC3);
    edge_dst = Mat::zeros(src.size(), CV_8UC1);
    this->TL = 50;
    this->TH = 150;
}

////////每次卷积均会丢失2*(k_size/2)|||此处为整除的边缘尺寸,
//void EdgeDetection::padding() {
//
//}

/**
 * @brief 判断一个点的坐标是否在图像内
 * @param r 行数
 * @param c 列数
 * @param rows 图像行数
 * @param cols 图像列数
 * @return
 */
bool checkInRang(int r, int c, int rows, int cols) {
    if (r >= 0 && r < rows && c >= 0 && c < cols)
        return true;
    else
        return false;
}

/**
 * @brief 从确定边缘点出发，延长边缘
 * @param maxEdgeMag 最大梯度边缘
 * @param edge 边缘结果图
 * @param TL 低阈值
 * @param r 行序号
 * @param c 列序号
 * @param rows 行数
 * @param cols 列数
 */
void trace(const Mat &maxEdgeMag, Mat &edge, float TL, int r, int c, int rows, int cols) {
    if (edge.ptr<uchar>(r)[c] == 0) {
        edge.ptr<uchar>(r)[c] = 255;
        // 遍历周围八个邻点
        for (int i = -1; i <= 1; ++i) {
            for (int j = -1; j <= 1; ++j) {
                float mag = maxEdgeMag.ptr<float>(r + i)[c + j];
                if (checkInRang(r + i, c + j, rows, cols) && mag >= TL)
                    trace(maxEdgeMag, edge, TL, r + i, c + j, rows, cols);
            }
        }
    }
}

/**
 * @brief 计算梯度方向以及非极大值抑制
 * @param sobel sobel算子对象
 * @param maxEdgeMag 最大梯度边缘图像
 */
void EdgeDetection::NonMaximumSuppression(const _Sobel &sobel, Mat &maxEdgeMag) const {
    for (int r = 1; r < src.rows - 1; ++r) {
        for (int c = 1; c < src.cols - 1; ++c) {/*避开边缘防止访问邻域时越界*/
            float x = sobel.dst_X.ptr<float>(r)[c];
            float y = sobel.dst_Y.ptr<float>(r)[c];
            float angle = atan2f(y, x) / CV_PI * 180; //当前位置梯度方向
            uchar mag = sobel.edgeMag.ptr<uchar>(r)[c];  //当前位置梯度幅值
            //垂直边缘--梯度方向为水平方向-3*3邻域内左右方向比较
            if (abs(angle) < 22.5 || abs(angle) > 157.5) {
                uchar left = sobel.edgeMag.ptr<uchar>(r)[c - 1];
                uchar right = sobel.edgeMag.ptr<uchar>(r)[c + 1];
                if (mag >= left && mag >= right)
                    maxEdgeMag.ptr<uchar>(r)[c] = mag;
            }

            //水平边缘--梯度方向为垂直方向-3*3邻域内上下方向比较
            if ((angle >= 67.5 && angle <= 112.5) || (angle >= -112.5 && angle <= -67.5)) {
                uchar top = sobel.edgeMag.ptr<uchar>(r - 1)[c];
                uchar down = sobel.edgeMag.ptr<uchar>(r + 1)[c];
                if (mag >= top && mag >= down)
                    maxEdgeMag.ptr<uchar>(r)[c] = mag;
            }

            //+45°边缘--梯度方向为其正交方向-3*3邻域内右上左下方向比较
            if ((angle > 112.5 && angle <= 157.5) || (angle > -67.5 && angle <= -22.5)) {
                uchar right_top = sobel.edgeMag.ptr<uchar>(r - 1)[c + 1];
                uchar left_down = sobel.edgeMag.ptr<uchar>(r + 1)[c - 1];
                if (mag >= right_top && mag >= left_down)
                    maxEdgeMag.ptr<uchar>(r)[c] = mag;
            }
            //+135°边缘--梯度方向为其正交方向-3*3邻域内右下左上方向比较
            if ((angle >= 22.5 && angle < 67.5) || (angle >= -157.5 && angle < -112.5)) {
                uchar left_top = sobel.edgeMag.ptr<uchar>(r - 1)[c - 1];
                uchar right_down = sobel.edgeMag.ptr<uchar>(r + 1)[c + 1];
                if (mag >= left_top && mag >= right_down)
                    maxEdgeMag.ptr<uchar>(r)[c] = mag;
            }
        }
    }
}

/**
 * @brief 双阈值处理
 * @param maxEdgeMag 最大梯度边缘图像
 */
void EdgeDetection::DualThresholdProcessing(const Mat &maxEdgeMag) {
    for (int r = 1; r < src.rows - 1; ++r) {
        for (int c = 1; c < src.cols - 1; ++c) {
            uchar mag = maxEdgeMag.ptr<uchar>(r)[c];
            //大于高阈值，为确定边缘点
            if (mag >= TH)
                trace(maxEdgeMag, edge_dst, TL, r, c, src.rows, src.cols);
                //小于低阈值，直接略过
            else if (mag < TL)
                edge_dst.ptr<uchar>(r)[c] = 0;
        }
    }
}

/**
 * @brief 边缘检测主程序
 */
void EdgeDetection::Edge_Canny() {
    cvtColor(src, src, COLOR_BGR2GRAY);//转换灰度图
    GaussianBlur(src, src, Size(5, 5), 0.8);    //高斯滤波
//    cv::Sobel(src, dst_X, CV_64F, 1, 0);//x方向灰度变化矩阵
//    cv::Sobel(src, dst_Y, CV_64F, 0, 1);//y方向灰度变化矩阵
    _Sobel sobel(src, k_size);
    //非极大值抑制
    Mat maxEdgeMag = Mat::zeros(src.size(), CV_8UC1);
    this->NonMaximumSuppression(sobel, maxEdgeMag);
    //双阈值处理及边缘连接
    this->DualThresholdProcessing(maxEdgeMag);
    imshow("dst", edge_dst);
    waitKey(0);
}
